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The Christodoulou memory is a nonlinear contribution to the gravitational-wave field that is 
sourced by the gravitational-wave stress-energy tensor. For quasicircular, inspiralling binaries, the 
Christodoulou memory produces a growing, nonoscillatory change in the gravitational-wave "plus" 
polarization, resulting in the permanent displacement of a pair of freely- falling test masses after the 
wave has passed. In addition to its nonoscillatory behavior, the Christodoulou memory is interesting 
because even though it originates from 2.5 post-Newtonian (PN) order multipole interactions, it 
affects the waveform at leading (Newtonian/quadrupole) order. The memory is also potentially 
detectable in binary black-hole mergers. While the oscillatory pieces of the gravitational-wave 
polarizations for quasicircular, inspiralling compact binaries have been computed to 3PN order, the 
memory contribution to the polarizations has only been calculated to leading order (the next-to- 
leading order 0.5PN term has previously been shown to vanish). Here the calculation of the memory 
for quasicircular, inspiralling binaries is extended to 3PN order. While the angular dependence of 
the memory remains qualitatively unchanged, the PN correction terms tend to reduce the memory's 
magnitude. Explicit expressions are given for the memory contributions to the plus polarization and 
the spin-weighted spherical-harmonic modes of the metric and curvature perturbations. Combined 
with the recent results of Blanchet et al. [Class. Quantum Grav. 25, 165003 (2008)], this completes 
the waveform polarizations to 3PN order. This paper also discusses: (i) the difficulties in extracting 
the memory from numerical relativity simulations, (ii) other nonoscillatory effects that enter the 
waveform polarizations at high PN orders, and (iii) issues concerning the observability of the memory 
in gravitational- wave detectors. 

PACS numbers: 04.25.Nx, 04.30.Db, 04.30.Tv, 95.30.Sf 



I. INTRODUCTION 

The primary purpose of this paper is to compute the 
post-Newtonian (PN) corrections to the nonlinear mem- 
ory piece of the gravitational-wave (GW) polarizations for 
quasicircular 1 , inspiralling binaries. We begin by review- 
ing the gravitational-wave memory (Sec. II A[) and then 
motivate its further study (Sec. II B[) . A summary of re- 
sults is given in Sec. II CI 



A. What is gravitational- wave memory? 

Gravitational-wave memory refers to the permanent 
displacement of an "ideal" GW detector after the GW 
has passed. An ideal detector is one which is only sen- 
sitive to gravitational forces — e.g., a ring of freely- falling 
test-masses — and is isolated from local tidal interactions. 
After the passage of a GW without memory, the detec- 
tor returns to its initial state of internal displacement (its 
state long before the passage of the wave). After the pas- 
sage of a GW with memory, the initial and final displace- 
ment states differ. 

There are two types of GW memory: linear and nonlin- 
ear. The linear memory has been known since the 1970s 
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1 quasicircular means that the binary orbit is circular up 2.5PN 
order radiation-reaction effects that cause the binary to slowly 
inspiral. 



[ESS] and arises from gravitational sources that produce 
a net change in the time derivatives of one or more of their 
source-multipole moments. For example, at leading order 
in a PN expansion the linear memory causes a net change 
in the GW field given by 



(i.i) 



where 2^ is the source mass-quadrupole moment, R is the 
distance to the source, TT means to take the transverse- 
traceless projection [Eq. (|2.2p ], and the A refers to the 
difference between late and early times. Geometric units 
(G = c = 1) are used here and throughout. 

A simple example of a source with linear memory is a 
binary on a hyperbolic orbit (gravitational two-body scat- 
tering) 0, H, @] . To see this consider the TT piece of the 
second-time derivative of the mass-quadrupole moment at 
Newtonian order and in the centcr-of-mass frame: 
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where fi — rrixm^jM is the reduced mass, M — mi + 
77i2 is the total mass, Xj is the relative orbital separation 
vector with magnitude r, Xj is the relative orbital velocity, 
and we have replaced second-time derivatives with the 
equation of motion Xj — —Mxj/r 3 . If we take t = to be 
the time of closest approach, then at very early (t — > — oo) 
and very late (t — > +00) times the relative velocities will 
be finite (tJco = y / 2E/fi where E is the orbital energy), 
while the second term in Eq. (|1.2p falls off like 1/r and 
is negligible. The memory for a hyperbolic binary thus 
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arises from the difference in the direction of the velocity 
vectors at late and early times: 
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Other systems with linear memory are those whose 
components change from being bound to unbound (or 
vice versa). These include binaries whose components are 
captured, disrupted, or undergo mass loss. Gravitational- 
waves with linear memory have been studied in the con- 
text of supernova explosions and their resulting neutron 
star kicks 0, H, H, [lO, ITU 113. lip, asymmetric mass loss 
due to neutrino emission [7L Il3l \l4. fla . llH ] (see Ref. fl7l 
for a recent review), or gamma-ray-burst jets [H, Hfll • 
Linear memory can also arise from GW recoil in binary 
black-hole mergers 21]. Sources whose components re- 
main bound generally do not display linear memory (but 
see Sec. IVBI for a caveat). 

A general formula for the linear memory produced by a 
system of N bodies with changing masses Ma or velocities 
v A is given by Thorne [H (also Ref. 0): 
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where the masses are unbound in their initial or final 
states (or both), the A means to take the difference be- 
tween the final and initial values of the summation, and 
N is a unit vector that points from the source to the ob- 
server. For example, the masses Ma and velocities va 
might refer to two or more particles on a scattering orbit 
(as discussed above), or they might refer to the various 
pieces of a star that becomes unbound. This formula is 
essentially the standard Lienard-Wiechert solution of the 
space-space part of the linearized Einstein equations (with 
a source term given by the stress-energy tensor of N non- 
interacting point particles). 

The nonlinear memory was discovered inde pen dently 
by Payne [23|, Blanchet and Damour [24], |25| . and 
Christodoulou [26j]. It is often referred to as the 
"Christodoulou memory." The nonlinear memory 2 arises 
from a change in the radiative-multipole moments that is 
sourced by the energy flux of the radiated gravitational 
waves. One can heuristically understand the origin of the 
nonlinear memory as follows: Consider the Einstein field 
equations in harmonic gauge [27| : 

□ ft"' 3 = -16n(-g)(T^ + tj£) - ft°^ft^ + h^h a ^ vl 

(1.5a) 



where 
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is the gravitational- field tensor 3 , g is the determinant of 
the metric g a 8, T a @ is the matter stress-energy tensor, 



f is the Landau-Lifshitz pseudotensor, □ = — d\ r + V is 
(1.3) the flat-space wave operator, a comma denotes a partial 



derivative ( ^ = <9 M ), and V is the flat-space Laplacian. 
While there arc many nonlinear terms that enter the right- 
hand side of Eq. (|1.5ap , there is a particular piece of f£ L 
[see the last term of Eq. (2.7) in Ref. [29(] that is equal to 
the GW stress-energy tensor [3(1 Hl| : 
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is the GW energy flux, 
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is a unit radial 



wnere lidn la w " ^^-bj 11 ttj j 
vector, the angle-brackets mean to average over several 

wavelengths, ftJ T = ft„ T , and we have used the plane- 
wave approximation hj^- as F a b(t—R)/R+0(R~ 2 ). When 
applying the standard Green's function to the right-hand 
side of Eq. (|1.5ap , this piece — proportional to the energy 
flux of the emitted GWs — yields the following correction 
term to the GW field HI: 
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where Tr is the retarded time. This equation shows that 
part of the distant GW field is sourced by the loss of GW 
energy. Note that while the magnitude of the nonlinear 
memory approximately scales with the total radiated GW 
energy Ai? gw , the angular-dependent unit vectors in the 
integrand of Eq. (|1.8p imply that the nonlinear memory is 
not directly proportional to AE gw . The memory should 
not be mistaken as a change in the monopolar piece of the 
1/R expansion of the metric. Rather it is a change in the 
quadrupolar (and higher-order) pieces of the 1 /R part of 
the TT projection of the asymptotic spatial metric (see 
Sec. ED] below). 

Since the nonlinear memory occurs in any system that 
radiates GWs, systems that are usually considered to have 
vanishing linear memory (such as bound binaries) have a 
nonvanishing nonlinear memory. Thorne [13] has shown 
that the nonlinear memory [Eq. (|1.8p ] can be described 
by his formula for the linear memory [Eq. (|1.4p ] if the un- 
bound objects in the system are taken to be the individual 
gravitons with energies Ea — Ma/(1 — v 2 ^) 1 / 2 and veloc- 
ities v J A = cn 3 A . The correspondence between these two 
equations also holds for null sources that contribute only 
to the linear memory. For example, the linear memory 
from a massless neutrino (or any other null particle) can 
be described either as a discrete sum over the memory 
from each individual particle [Eq. (|1.4p ], or by replacing 
the GW energy flux in Eq. (|1.8p with the energy flux of 
neutrinos 11511. 



2 Throughout this paper the terms "nonlinear memory" and 
"Christodoulou memory" are used interchangeably. When not 
otherwise specified, "memory" refers to the "Christodoulou mem- 
ory" . 



3 Different notations for h a P are used in the literature. For ex- 
ample, Blanchet [H and Will and Wiseman (WW) [H| both 
use the symbol h a P , but their conventions differ by a sign: 
h a P = Tlwtw = — ^b^* However, all agree on the sign conven- 
tion for the waveform hj^ . 
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The leading-order PN expansion of Eq. (|1.8[) is propor- 
tional to [see Eqs. (|2T]) . (|23T|l . and (pT22"|) ] 



(3) 

where X„- is the third-time derivative of the source- 

l j 

quadrupolc moment. If we specialize to quasicircular, 
inspiralling binaries with total mass M — mi + 7712, re- 
duced mass ratio r\ ~ mimi/M 2 , and angular orbital fre- 
quency uj(t), the integrand of Eq. JOJ) contains oscillatory 
terms proportional to r/ 2 (Moj) 10 / 3 e ±4lwt and nonoscilla- 
tory terms proportional to rj 2 (Muj) 10 / 3 . When performing 
the time integration, the oscillatory terms are effectively 
multiplied by the orbital timescale oc 1/cj: 
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while the nonoscillatory terms are effectively multiplied by 
the radiation-reaction timescale T rr oc (M/ri)(Mu))~ s / 3 : 
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This shows that while the oscillatory pieces of ShJ^ 
are a 2.5PN correction to the waveform amplitudes, the 
nonoscillatory piece enters at the same order as the famil- 
iar quadrupole-order piece of the waveform: 4 
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Motivation 
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The Christodoulou memory is a unique and interest- 
ing manifestation of the nonlinearity of general relativity. 
Analytic computations of gravitational radiation involve 
many types of nonlinearities, the origins of which are often 
obscured by the complex, iterative algorithms involved in 
solving the Einstein field equations. The Christodoulou 
memory, however, has a clear physical interpretation: It 
arises from the loss of GW energy from the system and the 
effect of this loss (through the GW stress-energy tensor) 
on the system's radiative mass-multipole moments. For 



This also follows from the approximate scaling of the nonlin- 
ear memory with the total radiated energy, 5hJ^ ~ A_E gw /R. 
Since the energy radiated during the inspiral is equal to the 
change in the orbital energy, AE gw Rj (r]M/2)(M/r), the non- 
linear memory has the same scaling as the quadrupolar waveform 
TT (0) 

hjf, oc (r/M/R)(M/r), where r is the orbital separation. 



comparison, GW tails are another interesting nonlinear, 
general-relativistic effect (see Sec. 3.4 of Ref. [33j for refer- 
ences). They arise from the last term in Eq. (|1.5ap . which 
modifies the flat-spacetime wave operator on the left-hand 
side and causes backscattering of the gravitational radia- 
tion as it propagates through the curved spacetime around 
the binary [34]]. Tail effects enter the waveform at 1.5PN 
order. The Christodoulou memory, on the other hand, 
arises from nonlinear interactions at 2.5PN order, but 
affects the gravitational waveform at leading (OPN) or- 
der. Both tails and the Christodoulou memory are hered- 
itary — their contribution to the GW field at any given 
retarded time depends on the entire past history of the 
source [see e.g., Eq. (j2.21[) ]. However, the Christodoulou 
memory is especially sensitive to the motion of the source 
in the distant past (see Sec. IV C|) , while tails are primarily 
sensitive to the recent past (see also Sec. 4 of Ref. |35j|). 
Tails and most other nonlinear effects that enter the wave- 
form cause oscillatory corrections to both the + and x 
polarizations. In contrast, the Christodoulou memory 
causes a nonoscillatory shift in the amplitude of the + 
polarization only 5 . This amplitude shift starts small at 
early times (when the binary is widely separated) and 
slowly grows during the inspiral. As the binary compo- 
nents merge the memory rapidly grows and then saturates 
to a final value during the ringdown phase. The details 
of how the me mory reaches its saturation value are ex- 
plored in Refs. 
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|36| |. Unlike tails, the nonoscillatory 
pieces of the memory do not affect the orbital phase of a 
quasicircular binary up to 3.5PN order (the highest order 
to which the phase has been computed), but they could 
modify the phase at higher PN orders 6 . 

Over the past three decades the PN corrections to 
the oscillatory piec es of the GW polarizations have been 
computed [H, l37l [H, EE] , most recently at the 3PN 
order [40| . These computations were motivated by 
the development — and finally the operation — of a global 
network of grou nd-based GW interferometric detectors 
[4ll l42l [43|, 13]. These detectors are primarily sensi- 
tive to the oscillatory components of the GW signal, as 
most of the signal power lies in the lowest-order oscilla- 
tory modes of the radiation (for nearly circular orbits). 
More recently, successes in numerical relativity (NR) 

EE m, S3, m IS EE El El H M HE HE M m, m 

have necessitated the need for accurate PN waveforms to 
compare with the results of binary black- hole (BH) merger 
simulations @, M, M, HI, M, HE ES HI M, H Ell, 



5 This statement is only true for standard choices of the polarization 
triad (see Sec. Ill All . Since a rotation of this triad by an angle 
about the propagation direction transforms the GW polarizations 



via hj 



[h+ — ih x )e , a purely + polarized wave can 
become mixed-polarized. 

For example (and using notation defined later), the GW luminos- 
ity C = —E oc y\ |^im| 2 w iH have a term proportional to |/i2o| 2 
which will affect the luminosity and hence the phase at relative 
5PN order (see Sec. lIII Bl below for scalings). However, one should 
note that the oscillatory terms that arise from the memory inte- 
gral [cf. Eq. Ill.l0a|l and associated discussion] affect the phasing 
beginning at relative 2.5PN order. 
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izaizaiziizaizflizaizaizai- 

Computations of the nonlinear memory's contribution 
to the waveform have not progressed as far. Wiseman and 
Will [32| first calculated the nonlinear memory's leading- 
order effect on the waveform polarizations for a quasicir- 
cular, inspiralling binary. Their result was confirmed in 
Refs. (35l. kfi l80l] . The 0.5PN corrections to the memory 
were calculated in Ref. [Z(| and found to vanish. The pri- 
mary purpose of this paper is to compute the corrections 
to the leading-order formula for the nonlinear memory to 
3PN order. 

Why are these corrections needed? First, while GW in- 
terferometers are mostly sensitive to the oscillatory parts 
of the GW, the memory piece of the signal could be de- 
tectable for certain sources. During a signal's observation 
time, the nonlinear memory causes a growing change in 
the signal's amplitude. This change contributes power at 
low frequencies. For many sources this power is swamped 
by the detector's low-frequency noise. But for sources 
with large signal-to-noise ratios — especially supermassive 
black- hole binaries in the low- frequency LISA [81] band — 
the memory can be detectable. For example, the memory 
from the merger of two 10 6 M Q black holes should be de- 
tectable by LISA out to a redshift of z » 2 (see Fig. 3 
of Ref. 36]). The detectability of the nonlinear memory 
was previously considered by Thorne [22| and Kennefick 
[80| and is discussed further in Sec. IV Dl below and in 
Refs. [Ill, [36|. Central to the issue of detecting the mem- 
ory is knowing its magnitude, and computing higher-PN 
corrections to the memory will help to determine this more 
accurately. 

Second, the comparison of PN waveforms with numer- 
ical relativity simulations could also benefit from more 
accurate expressions for the memory. While numeri- 
cal relativity accounts for all of the nonlinear effects of 
general relativity, it is difficult to compute the mem- 
ory accurately in these simulations: The modes of the 
waveform that have memory are very small and depend 
sensitively on the initial separation of the binary (see 
Sec. IV CI for details). The PN corrections to the mem- 
ory computed here could supply initial conditions for the 
waveform modes in the numerical simulations. As has 
been done with the oscillat ory pieces of the waveform 

IsH lE^, IBS, IBS, IEH> IHH> lsB> IBa, I^S, [^3, [zH Iz^, I^S, ) , it 

will also be insightful to compare the memory calculated 
in future simulations with the PN expressions computed 
here. 

Lastly, it is aesthetically pleasing to have all of the 
pieces — both oscillatory and nonoscillatory — of the wave- 
form amplitudes computed consistently to the same post- 
Newtonian order. Combined with the results of Ref. [io| . 
this work completes the waveform polarization amplitudes 
to 3PN order. 



7 The factor of 2 discrepancy in the memory waveform computed 
in Ref. l32l is either a typo or due to a different choice of nor- 
malization for the polarization tensors in Eq. i l'2.:il below. Note 
that the relative amplitudes shown in Fig. 1 of Ref. |32j have the 
correct values at large separation. 



C. Summary 

The bulk of this paper is devoted to calculating the 
nonlinear memory's contribution to the + waveform po- 
larization and its PN corrections. Readers uninterested in 
the details should skip to the results presented in Sec. lIVI 
Readers uninterested in lengthy formulas can skip directly 
to Figs. [1] and [21 the discussion (Sec. EJ, and conclusions 
(Sec. I VI I) . An outline of the paper and a summary of its 
main results are presented here: 

The necessary PN wave-generation formalism is briefly 
reviewed in Sec. HU Section ["II Al gives expressions for the 
waveform polarizations /i+, x m terms of the radiative- 
multipolc moments. Expressions for the spin-weighted 
spherical-harmonic mode decomposition of the waveform 
in terms of the "scalar" radiative-multipole moments 
are also introduced. Post-Newtonian calculations of the 
radiative-multipole moments are usually performed en- 
tirely in terms of symmetric-trace-free (STF) tensors, and 
are then (if needed) decomposed into spherical-harmonic 
modes at the end of the calculation. However, the com- 
putation of the memory is significantly easier if one works 
almost entirely from the start in terms of "scalar" quan- 
tities that are the coefficients of tensor quantities decom- 
posed on the basis of scalar or spin-weighted spherical 
harmonics. Sections III Bl and III CI review the multipolar- 
post-Minkowskian (MPM) formalism which relates the 
radiative-multipole moments Ml, and Vl to the canoni- 
cal moments Ml and Sl and the source moments Zj, and 
Jl- Section Hi Dl focuses on the hereditary contributions 
to the radiative multipoles: tails and memory. The re- 
sults of Blanchet and Damour 25j, which serve as the 
starting point of the memory calculation, are reviewed. 
The main result of this section is Eq. (|2.31|) . which ex- 
presses the memory piece of the radiative mass-multipole 
moments in terms of time and angular integrals over the 
GW energy flux from the source. 

The main part of the memory calculation — explicitly 
evaluating the nonlinear memory contributions to the 
radiative-mass multipoles for quasicircular binaries — is 
detailed in Sec. IIIII The computation of the angular in- 
tegrals is discussed in Sec. MI Al Section IIIIBI discusses 
which (l,m) modes need to be calculated to determine 
the memory contribution to the waveform to 3PN order. 
In Sec. IIII CI we discuss how the time integral over the 
past history of the source is computed. In particular we 
need to use a model for the adiabatic inspiral that is ac- 
curate to at least 3PN order. This model is contained 
in Eq. (|3.19|) . which gives the time evolution of the PN 
parameter x = (Mu;) 2 / 3 to 3.5PN order. This formula 
(which easily follows from the results of Refs. [H, |83| |) 
expresses the well-known result for the 3.5PN frequency 
evolution in a form that is useful for computations involv- 
ing the PN parameter x. 

The main results of this work are listed in Sec. IIVI the 
3PN memory contributions to the + polarization h+ of the 
gravitational waveform and its spin-weighted spherical- 
harmonic modes hi m . Those expressions can be directly 
combined with the oscillatory pieces of the waveform given 
to 3PN order in Ref. [4fJ . Figure Q] illustrates the mem- 
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ory's angular dependence and the relative sizes of the 
various PN corrections. The PN corrections have little 
effect on the memory's angular dependence, but tend to 
decrease the overall magnitude of the memory. 

Aside from the Christodoulou memory there are ad- 
ditional nonoscillatory (DC or "direct current") time- 
varying terms present in the waveform. Section IV Al dis- 
cusses a new type of nonlinear, nonhereditary DC term 
discovered by Arun et al. (35|. Unlike the Christodoulou 
memory, this nonlinear, nonhereditary DC term affects 
the x polarization and first enters the waveform at 2.5PN 
order. Section IV Bl discusses an additional class of lin- 
ear, nonhereditary DC terms that enter the waveform at 
5PN and higher orders. These terms arise from the ef- 
fects of radiation-reaction or nongravitational forces on 
the source-multipole moments. 

Section IV CI discusses the challenges in extracting the 
nonlinear memory from numerical simulations of binary 
black holes. It will be difficult for current simulations 
to accurately extract the memory waveform. One of 
the reasons for this is illustrated in Fig. [21 Numerical 
relativity simulations can most accurately compute the 
(l,m) — (2,2) mode of the waveform; but the nonlin- 
ear memory is present only in the m = modes (for bi- 
naries orbiting in the x-y plane). For simulations that 
can directly compute the metric-perturbation modes hi m , 
the largest memory mode /120 is an order-of-magnitude 
smaller than the /122 mode. For the majority of simula- 
tions that compute the spin-weighted spherical-harmonic 
modes of the ^4 Weyl scalar, ipi m = hi m , the situation 
is significantly worse: The largest memory mode "020 is 
nearly three orders of magnitude smaller than the domi- 
nant V>22 mode in the late inspiral. There are also difficul- 
ties in determining the two integration constants needed 
when computing the hi m modes from the ipi m modes. The 
memory's strong dependence on the past history of the 
source also introduces significant errors unless the simula- 
tions start with very large binary separations. The results 
of this paper could help to alleviate some of these difficul- 
ties by providing initial conditions for the m = modes. 



Section IV Dl discusses the detection of the memory. 
While the memory from a signal in the distant past is 
unobservable, it is possible to detect the buildup of the 
memory from a passing GW. Previous work by Thorne 
[22] ] and Kennefick [8(| has estimated the nonlinear mem- 
ory's signal-to-noise ratio in laser interferometers. How- 
ever, the signal-to-noise ratio is sensitive to the details 
of how the memory rises to its saturation value, which 
previous works have not properly modeled. 

Conclusions and suggestions for further work are pre- 
sented in Sec. IVII Some results relegated to the appen- 
dices are potentially useful in other applications: Ap- 
pendix|X]gives an explicit formula for the angular integral 
of the product of three spin-weighted spherical harmonics. 
Appendix [B] gives an explicit prescription for calculating 
the Weyl scalar modes ipi m from the known expressions 
for him presented here and in Ref. (4p| . 



D. Notation 

This paper borrows notations and conventions from 
Kidder Thorne [27[, and Blanchet, et al. 0. We 
generally set G = c = 1, except in certain equations 
where we wish to make explicit the post-Newtonian or 
post-Minkowskian order of the various terms. We use the 
notation 0{n) to denote post-Newtonian (PN) correction 
terms of order 0(c~ n ). Spacetime indices are denoted 
with Greek letters, spatial indices with Latin letters. In 
Sec. Ill Bl we use Blanchet's [28[ notation for the gravita- 
tional field, h a/3 = —h a P (see footnote[3]). Multi-index no- 
tation is denoted with a capital subscript: Al = A^...^,. 
A multi-index L on a vector denotes a product of I vec- 
tors: xl — Xi 1 Xi 2 . . . Xi r Repeated spatial indices and 
multi-indices are summed regardless of their relative po- 
sitions. Symmetric-trace- free (STF) spatial tensors are 
denoted with capital script letters (as in Ul)- The corre- 
sponding "scalar" versions of these moments [their coef- 
ficients on the basis of the STF spherical-harmonic ten- 
sors y l £ n (defined below)] are denoted with nonscript cap- 
ital letters (as in Ui rn ). Spherical-harmonic indices {l,m) 
are raised or lowered arbitrarily (i.e., Ui m = U Lm ). Sym- 
metrization and STF projection are denoted by enclosing 
the relevant indices by () or <>, respectively. Indices 
that are left out of the symmetrization or STF projec- 
tion are displayed with an underbar: A^Bajy. Time 
derivatives are denoted by an overdot or by superscript 
parenthesis: M.^ = d p M. L / dT p . To avoid confusion 
with the azimuthal harmonic index m, we denote the sum 
of the binary's component masses as Me m\ + 7712. This 
is distinct from the total binary (ADM) mass or mass- 
monopole moment, which we denote M. Euler's constant 
is denoted by = 0.577 21 . . .. The symbol <f) denotes an 
azimuthal angle while <p denotes an orbital phase. 



II. POST-NEWTONIAN WAVE GENERATION 
FORMALISM 

A. Waveform polarizations and mode 
decomposition 

Following Kidder we begin by introducing an 

asymptotically-flat radiative coordinate system described 
by spherical coordinates (T, R, 9,$) and their orthonor- 
mal basis vectors (eT, e^, ee, e$). The coordinate origin 
is the center-of-mass of the source. The retarded time in 
radiative coordinates is Tr = T — R. 

The gravitational waveforms are usually expressed as 
the transverse-traceless (TT) piece of the asymptotic met- 
ric perturbation decomposed into a sum over radiative 
mass- and current-multipole moments: 



j tt _ 
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Here Tlijki is the TT projection operator 



n 



ijkl 



PikPjl 



1 



:PijPkl, 



(2.2) 



where = 5ij — NiNj. The radiative mass- and 
current-multipole moments Ml(Tr) and Vl(Tr) are STF- 
tensors, and, as discussed below, are related to the source- 
multipole moments which are defined in the near zone as 
integrals over the source. Equation (|2.1[) represents the 
most general outgoing- wave, transverse trace- free solution 
of the vacuum- wave equation □/i^- T = [27j . 

The plus and cross GW polarizations are related to the 
TT piece of the asymptotic metric via /£_(_ x — l^ij , 
where the TT-polarization tensors are defined in terms of 
a chosen orthonormal triad (N, P, Q): 



1 



^(PiPj QiQj): 

sZ = ^(P i Qj + P j Qi). 



(2.3a) 
(2.3b) 



The unit vector N points from the source to the observer. 
While in some calculations it is convenient to choose this 
direction to lie along the z axis ez, here it points in a 
general direction specified by the spherical polar angles 
(0,$). We choose the remaining vectors in accordance 



with Kidder's [84j notation: P = e@ and Q 
that our choice of normalization implies 



e$° Note 



e A e B 



for A, B = + or x , and 



hf? = 2(h + e\ 



(2.4) 



(2.5) 



While waveform calculations usually involve computing 
the STF radiative multipoles Ul and Vl and then using 
Eqs. ()2.ip - ()2.3p to compute the + and x polarizations, 
it can be more convenient to compute the polarizations 
by performing a mode decomposition of the combination 
h+ — ih x via 



h + — ih x = 



where 



■Am 



G 



EE< 

1=2 m=-l 



U lm {T R ) - -V lm (T R ) 



_ 2 r im (e,$) 



V2Rc 



.1+2 



Ira i 



(2.6) 



(2.7) 



Here the "scalar" mass and current multipoles are related 
to their STF counterparts by [Thorne [23], Eq. (4.7)] 



U 



! iii 



AiU L y l i 



Irn * 



(2.8a) 



This choice differs from that chosen by Blanchet and collabora- 
tors [35l. |40H . The Kidder [§31 and Blanchet 40] polarizations are 
related by x = — ft¥ x (© = t, * = t/2). There is also a dif- 
ference in the overall sign of hi m between the two conventions. 
However the definitions of the radiativc-multipole moments and 
hT7~ are the same in both conventions. 



V 



Im 



Bi v L y l [ 



lm * 



where 



16tt 



Br 



(2l+l)\\ 
32nl 



'(^ + + 2) 
21(1-1) ' 



a- 



(2Z + 1)!! V 2J(J + 1)(J-1)' 



(2.8b) 
(2.8c) 
(2.8d) 



and * denotes complex conjugation. 

The complex conjugates of these moments satisfy 
[Eq. (4.5) of Ref. 0]: 



jjlm* / -^ymjjl—m -y-lm* ^ -|\my£ — m 

The inverse relations of Eqs. (|2.8|) are given by 




2^ 1) \ ~» jjlm^jlm 

+ 1)0 + 2) ^ ' 



m= — l 



+ 1)0 + 2) ^. y ' 



(2.9) 



(2.10a) 



(2.10b) 



m=—L 



These relationships also hold for the other mass-type and 
current-type multipole moments that are used elsewhere 
in this paper. The y l ™ are the STF spherical harmonics 
and arc related to the familiar scalar spherical harmonics 
by 



i = n L = n <L>, 



(2.11) 



where n.j is a general unit radial vector. They are given 
explicitly by Eq. (2.12) of Thorne [13] or Eq. (39) of Iyer 
|85j ]. We also note that for nonspinning binaries the h lm 
satisfy [Eq. (78) of Ref. 0] 



h 



1-1 



l iln 



(2.12) 



The spin-weighted spherical harmonics are defined in 
terms of the Wigner d functions by 



,y /m (e,$) = {-iy 



21 + 1 

47T 



d l ms (Q)e l 



(2.13) 



Here 



d l ms (Q) = y/(l + m)\(l-m)Kl + S )\(l- S )\ 
X (-l) fe (sin f ) 2fc+s - m (cos |)«+m-.-afc 



k=k, 



k\(l + m-k)\(l 



+ k)\ 



(2.14) 



where fej = max(0,m — s) and kf — min(i + m, I — s). 
The complex conjugates of the spin-weighted spherical 
harmonics satisfy 



s Y lm * = (-1) 



s+rn 



l—i 



(2.15) 



Working with the "scalar" multipole moments U lm and 
V lm can be more convenient because it allows us to work 
with 21 + 1 scalars instead of the 21 + 1 independent com- 
ponents of an Z-index STF tensor. Using the "scalar" 
moments also allows angular integrals over the products 
of unit vectors to be expressed in terms of products of 
spin- weighted spherical harmonics, which are more easily 
evaluated using computer algebra programs. 
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B. The multipolar post-Minkowskian formalism I: 
relating the radiative and canonical moments 

One of the purposes of a gravitational wave-generation 
formalism is to relate the radiative-multipole moments 
IAl and Vl that appear in the wavezone expansion 
of hj? [Eq. (|2.ip ] to some other family of multipole 
moments — the scmrce-multipole moments — which are de- 
fined in terms of integrals over the stress-energy pseu- 
dotensor of the matter and gravitational fields of the 
source. One procedure for relating these families of 
multipole moments is the multipolar-post-Minkowskian 
(MPM) iteration scheme developed by Blanchet, Damour, 
Iyer, and collaborators. This method is briefly sum- 
marized here and reviewed in detail by Blanchet [28| 
(see also Ref. [33| for a dated but much shorter review). 
Other PN wave-generation formalisms are discussed in 
Refs. [H, HI, H3] - 

The first step in the MPM procedure is a post- 
Minkowskian (weak field) iteration of the Einstein field 
equations subject to the harmonic gauge condition 
[Eqs. Ql.Sp ], This involves expanding the metric deviation 
h al3 = \/—gg a P — ri a P — —h al3 in powers of the gravitation 
constant G: 



Ghf 



2, a/3 



G z h 



G n K 13 



(2.16) 



and substituting into the vacuum Einstein equations and 
the harmonic gauge condition, resulting in a system of 
wave equations, 



Uhf = Af[h u ... dph 



a/3 







(2.17) 



Here A"' 3 represents the appropriate expansion of the 
right-hand side of Eq. (fL5a|) with T afi = 0. 

The next step consists of performing a multipolar ex- 
pansion of the KfP — an expansion in L/r where L < r is 
the size of the source and r is the field point. The coef- 
ficients of the powers of L/r can be expressed in terms 
of a new family of multipole moments. At linear order 
Thorne [27J has shown that the most general solution for 
(valid outside the source and up to an infinitesimal 
gauge transformation that preserves the harmonic gauge 
condition) is given by a multipole expansion that depends 
on only two types of STF moments, Ml and Sl- These 
moments are referred to as the canonical or algorithmic 
mass- and current-multipole moments. They represent an 
intermediate family of moments in between the radiative 
and source moments. 

Starting with this linear solution, each higher-order so- 
lution is generated by substituting the lower-order pieces 
into the right-hand side of Eq. (|2.17p and solving the re- 
sulting wave equation. This results in solutions at each 
order that take the form of a multipole expansion depend- 
ing on the canonical moments, h^f = h^\M.j,^ 

Because of the singularity at r = in the multipole ex- 
pansion, the ordinary retarded Green's function operator 
□ ~J yields divergent integrals. Instead one has to use a 
regularization procedure that consists of multiplying h^f 
by a factor (r/ro) B , where ro is an arbitrary constant 
length scale and B is a complex number. Applying the 



□ ~ t operator, and taking the finite part of the Laurent 
series expansion about B — > yields the following solu- 
tion to the wave equation at each order n (see Ref. [28[ 
for details): 



,a/3 



FP □: 

B=0 



A a/3 



(2.18) 



In order to satisfy the harmonic gauge condition as well 
as the wave equation, an additional piece v"@, which is 
constructed from the divergence of u^f and is a homo- 
geneous solution of the wave equation, must be added to 
Eq. (|2~T51) to yield the full solution to Eq. ([2~T7D at order 
n: 



,a/3 



(2.19) 



This result is then transformed from harmonic coordinates 
(ci, x l ) to radiative coordinates (cT, X 1 ) via the transfor- 
mation 



T R = t-~ 



2GM , / r 
■ In I — 



0(G 2 ), 



X 1 = x l + 0(G 2 ), 



(2.20a) 
(2.20b) 



where ro is another arbitrary length scale which is usually 
taken to have the same value as the ro in Eq. (|2.18[) , and 
A4 is the mass-monopole moment 9 . 

Taking the TT piece of the result and comparing with 
Eq. (|2.ip allows one to read off the relations between the 
radiative- and canonical-multipole moments. The result 
of this procedure is listed (to the highest PN order yet 
completed) in Eqs. (5.4)-(5.8) of Ref. [40| |. For illustration 
we show here the result for the radiative mass quadrupole: 



9 For quasicircular binaries this is related to the sum of the point- 
particle masses M by Eq. (1) of Ref. [8g|| : 



M 



M[l-^x+^x 2 (9 + r 1 ) + 0(c- e ) 



where we have used Eq. (6.6) of Ref. 1401 to express the result 
in terms of the PN parameter x = (Ma;) 2 / 3 for orbital angular 
frequency u>. 
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U ij {T R ) = Mf){T R ) 



2GM 





h( 2r ) 


11" 






1 — oo 




+ 12 



y 7 c 5 



X 3 <^K 3 >>)dr 



G 
7s 



,GW 



In 2 



Tr-t 
2t 



70 



2r 



124 627" 



44100 



M f j > {T)dT + 0[j), (2.21) 



where the constant To = ro/ c. 

Let us examine the various types of terms that ap- 
pear in Eq. (|2.21|) : First is the leading-order instanta- 
neous term that depends directly on the retarded con- 
figuration of the source. This term is familiar from the 
standard quadrupole formalism and enters at leading or- 
der in G [although, as shown below, this term contains 
0(G) corrections when the A4ij are expressed in terms of 
the source moments]. The remaining terms on the first 
and second lines are the nonlinear correction terms that 
enter at 0(G 2 ) in the MPM iteration scheme. The first 
of these is the leading-order tail term, which affects the 
waveform at 1.5PN and higher orders. 10 Tails arise from 
the scattering of GWs off of the monopole moment A4 of 
the source [25j. The next term is the leading-order non- 
linear memory, which is the primary focus of this paper. 
Because of the integral over the infinite-past history of 
the source, both the tail and memory are called hered- 
itary terms [25j . On the second line we have nonlinear 
instantaneous terms that affect the waveform at 2.5PN 
and higher orders [8f|. The last term is the cubically- 
nonlincar tail-of-tails term which affects the waveform at 
3PN and higher orders [88j, l9fj • The other radiative mass- 
and current-multipole moments have analogous correction 
terms [40]. For the rest of this paper we will primarily 
concern ourselves with the terms on the first line, and, in 
particular, the nonlinear memory term and its corrections 
(see Sec. HIS] below). 



C. The multipolar post-Minkowskian formalism II: 
relating the canonical and source moments 

The canonical moments M l and Sl do not have sim- 
ple closed-form expressions in terms of integrals over the 
source [28|. Instead they serve as intermediate moments 
that are related to a family of six types of scwrce-multipole 
moments: {2L,JL,yVL,X L ,y L ,Z L }. The mass- and 
current-source moments Xl and Jl tend to dominate 
over the four remaining moments which enter as 2.5PN 
corrections. Expressions for all of the source moments 
as explicit integrals over the stress-energy pseudotensor 



Note however that the 11/12 term in the integral is actually a 
nonlinear instantaneous term proportional to MM.^ (Tr). The 
same is also true for the analogous term in the tail-of-tails integral. 



[t q/3 = (167r) 1 x the terms on the right-hand side of 
Eq. (Tj~5aj) ] are found in Eqs. (85)-(90) of Ref. [H. 

The procedure for relating the canonical and source mo- 
ments is discussed in Ref. [401 ] . To summarize, it consists 
of performing the same MPM iteration discussed above, 
except the multipole expansion of the metric deviation is 
in terms of the source moments instead of the canonical 
moments. The result is a solution h"^[lL, ■ ■ -Zl] at each 
post-Minkowskian order that is related to the canonical 
moment solution h^f[ML, Sl] by a gauge transformation. 
This relationship between the two metrics can be trans- 
lated into a relationship between the canonical and source 
moments: 

Ml =Tl + GSI l + 0(G 2 ), (2.22a) 
S L = Jl + G5J l + 0(G 2 ), (2.22b) 

where the correction terms 62 l and 8Jl are functions 
of the six source moments. These corrections modify the 
leading-order mass- and current-source moments at 2.5PN 
and higher orders and are given (up to 3PN order) by 
Eqs. (5.9)-(5.11) of Ref. [Z(J. 

The final step in the gravitational wave-generation for- 
malism consists of matching the wave-zone MPM expan- 
sion of the metric in terms of the source moments with a 
post-Newtonian near-zone solution of the nonvacuum Ein- 
stein equations. This matching takes place in the region 
outside but close to the source where both approximation 
schemes are valid and yields an explicit relationship be- 
tween the source moments and the PN expansion of the 
near-zone metric. Solving the PN equations of motion for 
the source and substituting back into the PN-expanded 
stress-energy pseudotensor T a/3 (which depends on the 
metric and matter stress-energy tensor) yields explicit ex- 
pressions for the source moments in terms of variables 
describing the source (see Ref. [28| and references therein 
for the many nontrivial details of this procedure) . When 
specialized to quasicircular inspiralling compact binaries, 
the resulting expressions for the source moments are listed 
up to 3PN order in Eqs. (5.12)-(5.25) of Ref. 0. 

D. Hereditary contributions to the 
radiative-multipole moments 

Blanchet and Damour [25| give general expressions for 
the leading post-Minkowskian order hereditary contribu- 
tions (tail + memory) to the radiative mass- and current- 
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multipole moments for arbitrary I: 



U L =M ( P 



GU 



(tail) 



GU 



(mem) 



0(G 2 



GV 



(tail) 



+ 0(G 2 ) + 0(G/c 5 ). 



0(G/c 5 ), 
(2.23a) 

(2.23b) 



The neglected terms at 0(G 2 ) are tail-of-tails like terms 
and other cubically nonlinear interactions. The neglected 
0(G/c 5 ) terms are quadrupole-quadrupole nonlinearities 
[as in the second line of Eq. (|2.21[) ] and other instan- 
taneous "canonical moment x canonical moment" type 
terms. 

Using updated notation and incorporating the leading- 
order synchronous terms (those depending directly on 
events on the past null cone) into the tail integrals, the 
tail terms are given by [see Eq. (98) of Ref. [28j|]: 



U 



(tail) 



V 



(tail) 



2GM 



2GM 



T n r 



hi 



T n r 



In 



Tr-t 
2t 

Tr-t 
2t 



717 



M 



(r) dr, 
(2.24a) 



(2.24b) 



where the constants Ki and wi are 



2l 2 + 51 + 4 



1(1 + !)(/ + 2) ' j^fc 



1-2 x 



7T; 



Z-l 
1(1 + 1) 



i-i 1 

E 

fe=i 
(2.25) 

The memory term only affects the mass moments at 0(G) 
and is given by [Eq. (2.43c) of Ref. 0]: 



U 



mcm) _ 2c'- 2 (2Z + l) 



Tn 



dt 



dE 
dtdil 



n<L> 



Here 



(2.26) 

dtdn is the GW energy flux, i is a dummy variable for 
the radiative coordinate time; rii is a general unit vector 
centered at the source that points in the direction of the 
spherical polar angles (9, <fi) [distinct from the direction 
Ni from the source to the observer and the corresponding 



angles (6,$) appearing in the polarization waveforms]; 
and the angular integral is over the angles (9, (f>). Equation 
(|2.26p is the primary starting point for computing the 
nonlinear memory. 

Since Ref. [25| 's derivation of Eq. (]2 . 26[) was completed 
to second-post-Minkowskian order, their expression for 
the energy flux in Eq. (|2.26p formally depended only on 
the metric at first-post-Minkowskian (linearized) order 
[see Eq. (2.11b) of Ref. [25] and associated derivation]. 
However it is clear from other derivations (26l . [32j that 
the memory depends on the full GW energy flux. So 
Ref. (25|'s equation for the memory [their Eq. (2.43c)] is 
naturally extended to higher post-Minkowskian orders by 
using the full energy flux to the highest-PN-order known. 

Note that we are ignoring the "linear memory" con- 
tribution to £/j™ em ) arising from changes in the deriva- 
tives of the canonical-mass moments M.j} . We also note 
that changes in the derivatives of the canonical-current 
moments <S^ lead to linear memory contributions to the 
radiative-current moments Vl [25]. The linear memory 
does not directly contribute to the Christodoulou mem- 
ory and vanishes for quasicircular, inspiralling binaries 
that remain bound in the infinite past (but see Sec. IV Bl 
below). For bound astrophysical binaries whose compo- 
nents were formed, captured, exchanged, or underwent 
mass loss long before the GW driven regime, the linear 
memory is negligible. Note also that while there is no non- 
linear, hereditary memory contribution to the radiative- 
current multipoles, there is a nonlinear, "nonhereditary" 
DC (nonoscillatory) effect that arises from the 1.5PN cor- 
rection to Vijk', see the discussion in Sec. IV Al below. 

The GW energy flux can be computed from the GW 
stress-energy tensor and is given by [271 ] 



dE, 



gw 



dtdfl 



R 2 T, 



I Ml 



32tt 



16tt 



(2.27) 

where the angled brackets mean to average over several 
wavelengths, and we have used Eq. (|2.5[) to arrive at the 
last equality. Using Eq. (|2 .6[) we can write the energy flux 
in terms of the hi m modes: 



J 



dE s 



dtdVL 



16tt 



EE E 

Z'=2 l"=2 m'=- 



l" 

E 

/' m" — — V 



(hi' m /h*„ m „)- 



2 Y lm (9,0)_ 2 Y 



(2.28) 



Alternatively we can substitute Eqs. (|2.1[) or (|2.7[) into Eqs. (|2.27[) or (|2.28p to give the energy flux in terms of the 
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radiative-multipole moments [Eq. (4.14) of Thorne 

A 1 ) 7 A 1 ) „..„... A7/W 7/W 



(/' + !)!(/" + !)! 



8/'' 



>!(;" + i)! afcp p 

or its simpler expression in terms of the "scalar" moments 



Ual-l V bL»-l n L'-inL»-l + UacL>-2 V bl"-2 n L'-2n L »-2 



(2.29) 



dE, 



gw 



eft 32tt 



1 /rr(l) rr(l) T E2,l'm' rpE2,l"m", y{l) T/ (l) 

■^Ott 2-^i \ U l'rn' u l"m" J -jk ± jk ' v l'm' v I" m" 



(!) 77(1) rpE2,l'm' rpE2,l"m"j_ T/ (l) T/ (l) rpB2,l' m' rpB2,l" m"^ 07T (1) T/ (l) rpE2 X m j,B2 1 



l' ,V' : ra " 



jk 



1 jk ' ZU l'm' V l"m" 1 jk 



(2.30) 



In these expressions the "pure-spin" tensor spherical 



harmonics T^ k 2 ' lm and T^f" im are related to the spin- 
weighted spherical harmonics by Eqs. (2.38) of Thorne 
[27l | . and the summation limits are as in Eq. (|2.28p . 

The memory contributions to the waveform polariza- 
tions are conveniently computed by directly evaluating 
the hi m modes of Eq. (|2.6[) in terms of the memory con- 
tributions to the Uim multipolcs. This requires computing 
the "scalar" version of the memory piece of the radiative 
mass-multipole moment U^™ era \ Combining Eqs. (|2.26[> . 
dUHlD, and (|2~TTj) yields 



nB2,lr< 



u, 



(mem) 
lm 



32tt / (I - 2)! 



„2-l 



dt / dn 



dE e 



c— y 2(Z + 2)! J-vTJ dtdSi^ Ylm ^' 

(2.31) 

This is the primary equation that we need to evaluate 
in order to compute the nonlinear memory and its PN 
corrections. 

For completeness, the combined contributions to the 
waveform modes from all of the hereditary and syn- 
chronous terms discussed in this section are given by 



2GM 



G 



V2Ri 



,1+2 



+ 



In 



2GM 



2r 



M ('+2) 

1V1 lm 



(r) 



lm 



(r) 



dr 



«^ +1) (T fl 



-7T,S (Z+1) 



(Tr) 



+ GU^ cm) (T R ) + 0(G 2 ) + 0(G/c 5 )|. (2.32) 

11 To translate from Thome's [27J notation for the radiative- 
multipole moments to the notation used here and in Ref. [84| , the 
following replacements are made in Thome's formulas: ' ! '2^4 ; — * 
Ul, (l) S Al -* V L , (i) ^m -* U lm , and WS Jm -> V lm . In this paper 
the symbols Xj, and 7; m refer to the source mass-multipole mo- 
ments, while Sl and S; m refer to the canonical current-multipole 
moments. 



III. EVALUATING THE MEMORY 
CONTRIBUTION TO THE RADIATIVE-MASS 
MULTIPOLES 

A. Computing angular integrals 

Computing the nonlinear memory's contribution to the 
waveform primarily consists of evaluating the radiative- 
mass multipoles in Eqs. (|2.26[) or (|2.31[) . The most in- 
volved step is computing the angular integral that appears 
in these expressions. If one uses Eq. (|2.26p and works with 
the expansion of the energy flux in terms of STF tensors 
[Eq. (|2.29p ] one needs to compute a sequence of angular 
integrals of the form 



dfl n l , 



(3.1) 



which have well-known expressions in terms of products 
of Kronecker deltas [Eq. (2.3) of Ref. [13]]. However, 
for large I values this approach involves expressions with 
many indices and a large number of terms from the sums 
in Eq. (|2.29[) . While it is possible to use tensor algebra 
software to perform the manipulations, it is much easier 
to evaluate the angular integral in Eq. (|2.3ip using the 
expression for the energy flux in Eq. (|2.28|) . In this case 
all of the angular integrals have the form 



dfl -0Y1 



2-1 I'm' 



v* 
2 1 l"m" 



(3.2) 



Using Eq. (|2.13|) one can easily derive an exact expres- 
sion for this integral (Appendix [X| . The result involves 
a complicated sum that is tedious to evaluate by hand. 
Computer algebra programs can easily compute the an- 
gular integrals in Eq. (|3.2[) by cither evaluating the sum in 
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the resulting analytic solution [Eq. (|A3|) ]. or by comput- 
ing the various harmonics via Eq. (|2. 13|) and performing 
the integrals symbolically. 

If we define the time derivative of the memory 
mass-multipole moment in Eq. (|2.3ip . C^/™ = 

dU ^ sm) /d T R , then combining Eqs. (f23T|) . (l2~28l) . ((A3]), 
and (|2.15p gives 



time integral in Eq. (|2.3ip leads to two types of integrals: 



+ (higher order terms) for m/0, and (3.7a) 



U, 



(mem)(l) 



'2(1-2)1 



oo oo 



I" 



EE E Em: 



m+m 



(1 + 2)1 ^ ^ ^ v 

v ' l'=2l"=2m'=-l' m"=-l" 



x{hi. m >ht,, m ,,)(%,T™ m ,_ m „_ m> (3.3) 

where Gf\ S2 , s i, „ „ is an angular integral related to 
Eq. (|3.2p and is given in Appendix |A"1 



B. Selection rules and time derivatives of the 
memory multipole moments to 3PN order 

To compute the nonlinear memory's contribution to the 
waveform polarizations to 3PN order, we must expand the 
energy flux in Eqs. (|2~3Tj) or to 3PN order. For a 

quasicircular binary with relative orbital speed v, orbital 
separation r, and orbital period T, the hi m modes for 
to 7^ have the leading-order PN scaling 



h, 



lm 



lm 

R 



77 M 
~R~r~ 



2+1 



R 



(3.4) 



Since the leading I = 2 term is of order 0(v 5 ), our knowl- 
edge of the waveform to relative 3PN order implies that 
we must evaluate the sums in the energy flux or Eq. (|3.3p 
up to maximum I' and I" values of l' max = Z max = 8. Ex- 
pressions for the him modes to 3PN order and / < 8 are 
given in Eqs. (9.3)-(9.4) of Ref. 0. The straightforward 
procedure for computing the time derivatives hi m is dis- 
cussed in Appendix [B] and uses the results of Sec. IIII CI 
below. 

When evaluating , what values of I and to do we 

need to compute? To determine this we first note that 
the (f> integrals in Eq. I|3.2p . 



r 2 ~ 2t) nc 

{ - r l'l"lm'-m"-m ^ 



'd(f>, (3.5) 



imply the selection rule to = to' — to". Since the hi m (for 
to 7^ 0) scale like 



h lm ex IxV+We-*"*, 



R 



the product appearing in Eq. 



(3.6) 



is proportional to 



x" dt 



— dx 

x 



M 



64(n — 4) T) 
(higher order terms) for to = 0, (3.7b) 



where x = 64r]x 5 / '(5Af) [1 + 0(2)] [see also Eqs. (fTTT0|) and 
Ref. [1H|]. The to ^ terms yield oscillatory contribu- 
tions to the waveform polarizations that enter at higher 
PN orders than the nonoscillatory, to = terms. In a 
complete computation of the waveform, these oscillatory 
pieces contribute to the full- waveform beginning at 2.5PN 
order. Since we are interested only in the nonoscillatory 
memory effect, we will only focus on computing the m = 
terms in U^ cm \ However the procedure discussed here 
can also be used to compute the m/0 oscillatory terms. 

We also note here that the restriction of the memory 
to the to = modes is largely a consequence of choosing 
our coordinate system such that the binary's motion is 
confined to the x-y plane. This choice allows the hi m 
modes to be proportional to e~ lmip ^ . A rotation of our 
coordinate system (or equivalently a rotation of the orbital 
angular momentum) would mix the to modes and lead to 
nonoscillatory memory terms in the to 7^ modes as well. 

To determine the maximum value of / that is needed 
in l^Q mcm \ we first consider the angular integral in 
Eq. (|2.26p . where the energy flux is given by Eq. (|2.29p . 
The STF properties of the radiative multipoles, combined 
with the Kronecker deltas that result from angular in- 
tegrals over the products of rij, require that the max- 



imum I for which hi 



or U, 



(mem) 



/' + I 
'max 1 ' 



II 

max' 



will be nonzero is 



For example, the first term in the 



energy flux in Eq. (|2.29p contributes to the U\ 
terms of the form: 



via 



ut em){1) oc y T 



nLnL'TiL"dQ 



(3.8) 



The STF properties of anc ^ Ml then require that the 
index coupling that maximizes I be of the form 



^(mem)(l) ^ yU^W^l) 



(3.9) 



implying l max = l' max + ^ max - Using the 3PN waveform as 
input (V = I" < 8), this implies that the angular integrals 
in Eq. (|3.3p will vanish for I > 16. However we note from 
Eqs. and (GO} that the PN scaling of C//™° m)(1) is 



U, 



i)U) 



77V'+ r + 6 , 



(3.10) 



where if is the orbital phase, x = (Mw) 2 / 3 is the standard 
PN expansion parameter for circular orbits, and n > 5. 
Combining with the above-mentioned selection rule, the 



which for l max = 16 (/' = I" = 8) implies a 6PN-order cor- 
rection relative to the leading-order f/^m ~ 0(v 10 ) 
term. Since our 3PN waveform allows us to consistently 
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compute corrections only to relative 3PN order, this im- 
plies that the maximum I up to which we can consistently 

compute the £^/™ em ^^ is l m &x — 10. 

Before proceeding with the computation of the 

, we note that the right-hand side of Eq. (|3.3|) 
is a function of the full h[ m modes [Eq. (|2.32p j. which are 

themselves functions of the memory multipoles t^/™ em ^ 
that we are trying to compute. However, the time deriva- 
tives t/^™ 0111 ^ 1 ) f the memory modes provide a very small 
contribution to the energy flux and can be neglected in 
comparison to the higher-order m ^ modes. For exam- 
ple while the oscillatory mode product in Eq. (|3.3j) scales 
like 



,10 



+ higher order terms, (3-11) 



the nonoscillatory (memory) mode product scales like 



hi'nh 



l'O n l"Q 



,,20 



E?' 



+ higher order terms. 



(3.12) 



Even at their lowest order, the /ij'o^z"o mode pairs con- 
tribute a relative order 0(v 10 ) (5PN) correction to the 
right-hand side of Eq. (|3 . 3[) and can be safely ignored. 
The oscillatory modes effectively act as generators for the 
nonoscillatory (memory) modes. 



The explicit computation of the uj^ em ^^ 



modes is performed using the computer algebra program Maple. First 
the sums appearing in Eq. (|3.3p are expanded, and the angular integrals are evaluated via Eq. (|A3p . The resulting 



expressions consist of long sums over the various hi>mih*,, m ,i terms. Next the explicit PN expansions for the hi m in the 
quasicircular case [given in Eqs. (9.3)-(9.4) of Ref. [40(] are explicitly differentiated with respect to time. The details 
of this time differentiation are discussed in detail in Appendix [Bj The essence of the calculation involves expressing 
the amplitude and phase of hi m entirely in terms of the PN parameter x and then computing the time derivative via 
hi m [x(t)] — x(dhi m /dx). The formula for x is itself derived in the next section from the 3PN GW luminosity and is 
directly related to the rate-of-change of the orbital frequency (which is itself a crucial quantity needed for the detection 
of inspiralling binaries with LIGO and other detectors). Once these derivatives are computed and substituted into the 
expanded sum on the right-hand side of Eq. (|3.3[) . the result is series expanded to the PN order appropriate for the 



given (l,m) mode. 12 The resulting expressions for the £/;™ e are given by: 



U. 



(mcm)(l) 



2(1 



256 
~21 



rj 2 x 5 



nx 



5/2 



2435 



1219 
"288" 

23 



144 12 



^+4nx^+, 

174 213 949 439 
1816 214 400 
126 714 689 41 



4447872 



48 



793 14 023 



1782 
16tt 2 



6336 
856 



V - 



4201 



3 105 
4168 379 



1584 
{2 lE + In 16x) 
, 142 471 , 



123 552 



■if 



46 332 



0(7)}, (3.13a) 
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+ 0(7U, (3.13b) 
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(3.13c) 



12 I have also checked that applying a similar procedure to the an- 
gular integral of the energy flux in Eq. I|2.28|l yields the correct 
formula for the GW luminosity to 3PN order [Eq. 1 13, 181 1 below] . 
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U 



10 



525 221 



15 752 880 V 385 



777 ar 



_ 79 841 784 
~ 9 979199 V 



198 570 240 
9 979199 



-V 



172 307 520 
9 979199 



rf + 0(2) 



(3.13e) 



The odd-^ moments vanish up to the required PN order: 



U. 



:',{} 



u, 



(mem)(l) 
50 



U. 



(mem)(l) 
70 



u, 



(mem) (1) 
90 



(3.14) 

Note that both {/ 2 ( ™ m)(1) and L/j nom)(1) are of order 
0(x 5 ). This is a consequence of the i max = l' max + Z^ax 
selection rule: If the energy flux is expanded only to lead- 
ing order (l' max = l'J[ iax — 2), the selection rules allow all 
I < 4. 



Computing time integrals over the past history 
of the source 



The final step in computing the memory contribution 
to the radiative mass-multipole moments is the evalua- 
tion of the time integral over the entire past history of 
the source [cf. Equation (|2.3ip ]. While these time inte- 
grals generally come in the two types listed in Eqs. (|3.7p . 
the multipolc components that contribute to the memory 

[L^™ em ^l involve only integrals of the second type. Arun 
et al. [35J discuss in detail an adiabatic model for the evo- 
lution of the source that consists of an inspiralling binary 
described by the leading-order (2.5PN) radiation-reaction 
formulas. Since we are here computing 3PN corrections 
to the memory terms considered in Ref. [35J, we must 
extend our adiabatic model of the inspiral to three PN 
orders beyond the leading-order model considered there. 

It is important to note that this adiabatic model is an 
idealization representing a binary whose orbit has been 
decaying along quasicircular orbits from infinite separa- 
tion solely via gravitational radiation reaction. In reality 
a binary is formed or captured with some initial separa- 
tion and eccentricity. This formation or capture, as well 
as perturbations to its orbital elements (e.g., via mass loss 
or gravitational 3-body scattering) may have caused addi- 
tional linear or nonlinear memory contributions. The adi- 
abatic model used here ignores these issues and assumes 



a perfect quasicircular inspiral from infinite initial orbital 
separation. This is likely to be a good approximation if 
deviations from quasicircularity occurred in the very dis- 
tant past (long before the binary is observed). Future 
work will test this approximation by including the effects 
of the binary's eccentricity (which grows in the past) [9l[ . 

The main "trick" in computing the time integrals of 
Eqs. (|3.13p is to simply change variables from time t to 
the PN expansion parameter x: 



[x(t)] n dt = 



x(Tr) x n 



dx, 



(3.15) 



where we have used the fact that x — > in the infinite past. 
While computing the Newtonian-order memory requires 
only the leading-order contribution to x, calculation of 
the 3PN memory requires the 3PN corrections to x. This 
3PN formula for x (which is often expressed in terms of Co 
and u> in the literature) is the essence of our 3PN adiabatic 
model. 

Applying the chain rule we can express x in terms of 
the orbital energy and the GW luminosity L = —E, 



dx (-£) 



(3.16) 



dt dE/dx' 

The 3.5PN orbital energy is given by Eq. (5) of Ref. [1 
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and the 3.5PN GW luminosity is given by Eq. (2) of 
Ref. SI: 
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(3.18) 
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Note that the 3.5PN term in Eq. (l3TT7|l is zero and we have substituted the values A = -1987/3080 and = -11 831/9240 
for the ambiguity parameters [92j, |93j, |9J] . Substituting these expressions into Eq. (|3.16[) and expanding to 3.5PN order 
yields 
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Using Eq. (|3.19|) in Eq. (|3. 15(1 and expanding the integrands to the appropriate PN order, the time integrals of Eqs. (|3.13p 
are easily computed. 13 



IV. RESULTS: MEMORY CONTRIBUTIONS TO THE WAVEFORM MODES AND POLARIZATIONS 

Rather than listing the resulting expressions for uj^ e , we instead list the memory contributions to the spin- weighted 
spherical-harmonic modes of the polarization waveform [Eq. (|2.7j) ]. These quantities are simply related via 
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l l0 
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5 R 



(4.1) 



where we have followed the notation of Sec. 9 of Ref. [40] ■ The notational parameter a accounts for the two commonly 
used choices for the polarization triad (see discussion in Sec. Ill Aj ): 



(+1) for the Kidder [84] convention, 
(—1) for the Blanchet et al. [40] convention. 

The resulting polarization modes in terms of Hiq are: 
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13 Note that when computing integrals and PN series expansions 
here and elsewhere in this paper, the natural-log terms ln(16x) 
and ln(x/xo) [cf. Equation (IB 2 1 1 are not treated as constants. 
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(4.3e) 



Choosing a = —1, these expressions can be directly combined with Eqs. (9.4) of Ref. [40J to give all of the 3PN 
contributions to the spin- weighted spherical-harmonic modes /i; m . Note that all of the contributions to the m = 0, 
even-? modes arise solely from the Christodoulou memory piece of Ui m . There are no hereditary- memory contributions 
to the radiative current-multipole moments Vi m . However, there is a nonhereditary DC contribution to /130 at 2.5PN 
order arising from the nonlinear corrections to Vjj&. This is discussed in Sec. IV Al below. 

Combining the above modes with Eq. (|2.6| . we can explicitly compute the memory contributions to the + waveform 
polarization. Following the notation in Ref. [40j we factor the waveform as 
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The memory contributions to H^ 1 ^ 2 ^ are: 
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where cq = cosO and sq = sinO. Choosing a = —1 and = t, these expressions can be directly combined with 
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FIG. 1: (color online). Dependence of the post-Newtonian (PN) corrections to the Christodoulou memory on binary inclination 
and orbital separation. The plots in the top row show the memory contribution to H+ [Eq. ((43}] at each cumulative PN order 
as a function of Q (the polar angle to the observer; = points along the binary's orbital angular momentum) for x — 1/5 and 
r\ — (left) and r\ — 0.25 (right). The different curves represent terms up to the following cumulative PN orders: solid (black) 
OPN; long-dashed (red) 1PN; short-dash-dotted (green) 2PN; long-dash-dotted (magenta) 2.5PN; short-dashed (blue) 3PN. The 
bottom row plots xH+ [see Eqs. (|4.4[) and (|4.5|l ] and shows the memory's dependence on the PN parameter x (which equals 
M/r at Newtonian order, where r is the orbital separation in harmonic coordinates). The labeling scheme is the same as in the 
top row. The PN corrections do not qualitatively change the angular dependence, but tend to decrease the magnitude of the 
memory. Since the 2.5PN correction vanishes for r\ = 0.25, the 2.5PN curve is identical to the 2PN one and is not displayed in 
the right column plots. For r\ — 0.25 the 3PN curve is nearly coincident with the 2PN curve. 



Eqs. (8.9) of Ref. [4(| to yield the total + polarization. 
For reference, we also recall the leading-order, nonmemory 
contributions to the polarizations: 

ftCQ) = ^2a^-x{\ + 4) cos 2(tp - $), (4.7a) 
R 

h { ° ] = -Aa\-xce sin2(<p- $), (4.7b) 
R 

where — l and $ = it/ 2 in the conventions of Ref. [4(ij |. 

We note the following features of the memory wave- 
form: (i) As pointed out previously by Ref. |4(|, there 
is no memory term at 0.5PN order; we see here that the 
1.5PN contribution to the memory also vanishes, (ii) The 
2.5PN term vanishes for equal-mass binaries. Even in the 
77 — > limit, the 2.5PN term has the smallest magnitude of 



all the nonvanishing PN terms computed here, (iii) Aside 
from the nonhereditary DC term discussed below, all of 
the hereditary memory terms only affect the + polariza- 
tion. This arises from the choice of our polarization triad, 
the fact that the i//™° m ^ are all real, and the absence of 
hereditary-memory contributions in the radiative-current 
multipoles. A rotation of the polarization triad would 
cause nonlinear memory contributions to both polariza- 
tions, (iv) As expected for m = modes in a planar 
system, the memory contribution to the polarizations is 
independent of the angle $. (v) Lastly, note that the 
memory pieces of the waveform are entirely free of any 
arbitrary constants (such as ro) that arise in the MPM 
formalism (Sec. Ill Bl) . 

Figure [T] indicates several more features of the PN cor- 
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rections to the memory waveform: modulo some factors, 
these plots display the memory contribution to the + 
waveform polarization [Eqs. (|4.4|) -(|4.6 |) ] from OPN to 3PN 
orders (e.g., the 3PN curves contain all PN terms up to 
and including the 3PN terms). In the top row we see 
that the angular dependence is qualitatively similar at 
each PN order, with the memory's amplitude peaking at 
O = 7r/2. The bottom plots show how the PN correc- 
tions to the memory depend on the PN parameter x, with 
x ranging from 1/30 [r ~ 30A/) to 1/5 (approximately 
corresponding to the last stable orbit at harmonic coor- 
dinate radius r » 5M). 14 While the OPN piece of the 
memory increases linearly with x, the higher PN correc- 
tion terms tend to decrease the memory: for example, 
for r\ — 0.25 and x — 1/5, the 3PN memory is smaller 
than the OPN memory by a factor of w 0.74. Although 
the memory contribution to h + vanishes in the 77 — > 
limit, the plots have factored out the leading-order r\ de- 
pendence and allow us to compare the equal-mass and 
extreme-mass-ratio limits. Notice that as the PN order is 
increased, the memory waveforms seem to be converging 
rapidly in the equal-mass limit, but converge more slowly 
for small mass-ratios. This is consistent with the behav- 
ior of the oscillatory pieces of the waveform and the GW 
luminosity. 



V. DISCUSSION 

Now that we have completed our derivation of the 
higher-order PN corrections to the memory waveform, we 
can address several other memory-related issues. The first 
two subsections below discuss additional nonoscillatory 
contributions to the waveform. In the last two subsections 
we address the challenges associated with computing the 
nonlinear memory in numerical relativity (NR) simula- 
tions and observing memory in gravitational-wave (GW) 
detectors. 

Throughout this paper and elsewhere in the litera- 
ture, memory effects are often referred to as DC (for 
"direct-current" ) effects in the waveform. This terminol- 
ogy should be taken to mean simply that the memory is a 
nonoscillatory modulation of the waveform, in contrast to 
the AC or oscillatory waveform modulations. In partic- 
ular we note that DC waveform corrections are generally 
not constant offsets in the polarizations but can vary with 
time. For bound binaries this variation usually proceeds 
on a slow (radiation-reaction) time scale rather than an 
orbital time scale. For example, the GW modes for qua- 
sicircular binaries can be expanded in a Fourier series of 
the form 



+00 



h lm (t)= c'w 



ifi(t)t 



(5.1) 



14 Note that x = (M/r)[l + 0(2)], where the PN correction terms 
can be inferred from Eq. (7.21) of Ref. [95| . 



where £1 = 1/t J uj(t)dt and Tij^ is some mode amplitude 



function. In the absence of radiation-reaction £1 and Ti}' 
are constant, but they are slowly evolving functions of 
time when radiation-reaction is included. In such a sum 
the DC terms are those corresponding to n = 0. We 
note that this Fourier index n is generally distinct from 
the azimuthal angular index m, but, for the conventional 
choice in which the orbit lies in the x-y plane, n = m, and 
we can remove the sum in Eq. (|5.ip . 



A. The nonlinear, nonhereditary DC term 

In addition to the nonlinear, Christodoulou memory 
which has been the primary focus of this investigation, 
a new type of nonlinear, zero- frequency term has recently 
been discovered by Arun et al. 35]. Unlike the hereditary 
memory discussed here, this term is nonhereditary and has 
its origin in the (mass quadrupole) x (mass quadrupolc) 
and (current dipole) x (mass quadrupole) terms that arise 
in the 1.5PN corrections to the radiative current-octupole 
moment [Eq. (5.6b) of Ref. 0]: 
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where 



Si = Ji — rjMtiabXaVb + 0(c 2 ), 



(5.3) 



and the last three terms in Eq. (|5.2[) give rise to the non- 
hereditary DC effect. 

We can gain further insight into this effect by expand- 
ing Eq. (|5.2[) on the basis of STF spherical harmonics. Us- 
ing Eq. (|2.8b[) to compute V 3m from Vijk and Eq. (|2.10aj) 
to express M.ij in terms of M% m , the radiative current- 
octupole moment can be written as 
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where we have used Si — J Z S* and M 2 ±i = (for or- 
bits in the x-y plane). Since for quasicircular and planar 
orbits the canonical-mass moments are proportional to 
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one can easily see that the sinu- 

(B) 



soidal dependence on tp cancels in Vg above. 

Unlike the Christodoulou memory the nonlinear, non- 
hereditary DC effect originating from modifies the x 
polarization waveform at the 2.5PN order. In the notation 
of Sec. IIV1 this term provides the following contribution 
to the him modes [Eq. (9.4g) of Ref. [H], 
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and to the x polarization [Eq. (5.10) of Ref. [3 
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It is not clear if there is a simple physical explanation 
for this nonlinear, nonhcreditary DC term. Because of 
its high PN order, it is likely to be of much less obser- 
vational significance than the Christodoulou memory. It 
is also not clear if this term leads to a "true" memory in 
a GW detector, i.e., a displacement in the detector that 
persists after the GW has passed [96]. It is possible that 
this nonhereditary, DC waveform correction grows dur- 
ing the inspiral phase, but then decays to zero during the 
merger and ringdown, leaving no net memory. Nonethe- 
less, it is interesting that this type of nonlinearity [the 
canonical-moment coupling in Eq. (|5.2p ] can give rise to a 
DC effect in the waveform. It is possible that other multi- 
pole interactions of this type could contribute additional 
nonhereditary DC terms at higher PN orders. 



The m/0 modes have time derivatives proportional to 
Ilm x e - ™ 1 ^*). For unbound orbits the time derivatives 
7^2 that enter the waveform can display linear memory 
from differences in the phase angle (p. 15 But for bound 
orbits, these modes yield purely oscillatory contributions 
to the polarizations. 

To see a linear DC effect for bound orbits one has to 
examine the to = modes. For quasicircular orbits the 
time derivatives of the to 7^ moments are dominated 
by the derivatives d n / ' dt n (e'^) oc u n e- %mv , with the 
derivatives of r(t) contributing small, oscillatory correc- 
tions at 2.5PN and higher orders. But for the to = 
modes, the derivatives of r(t) dominate. For example, for 
I = 2, to = 0, the second derivative of the source mass- 
quadrupole moment is 



4l ] (xrj 3 Mx 6 [l + 0(2)}. 



(5.11) 



This term arises from the noncircularity of the orbit in- 
duced by radiation-reaction. Analogous formulas hold for 
higher I values and for the source-current moments. Terms 
of this form vanish in the distant past, are slowly grow- 
ing, and nonoscillatory. They contribute DC terms to 
the waveform at 5PN and higher orders (and are thus 
negligible for most sources). Furthermore, these DC ef- 
fects in bound binaries need not only result from gravita- 
tional radiation-reaction. Consider, for example, a com- 
pact body orbiting inside a dense accretion disk that sur- 
rounds a massive black hole. Viscous forces inside the 
disk will also cause an inward radial velocity r, producing 
a slowly changing, nonoscillatory contribution to ij^ and 

Jjq , and a resulting DC contribution to the waveform. 
One can imagine other gravitational and nongravitational 
forces that might produce a similar DC effect. 

As with the nonlinear, nonhcreditary DC term dis- 
cussed in Sec. IV Al it is not clear if these linear DC ef- 
fects give rise to a permanent memory. It is possible that 
these DC contributions to the waveform will vanish at late 
times. 



B. Linear DC effects in bound binaries 

Linear memory is generally considered to arise from per- 
manent differences between late and early times in the 



derivatives of the source moments AZ^ and /S.j'^ 1 [or 



AD 



equivalently A/^ and AJ^]. While unbound orbits de- 
rive their memory from changes in the derivatives of the 
source moments, bound orbits are generally thought to 
have vanishing linear memory. Here we argue that bound 
orbits that undergo gravitational radiation-reaction — or 
any other nonperiodic, secular change — can also display 
nonoscillatory (DC) waveform components resulting from 
long-timescale changes in the derivatives of the source- 
multipolc moments. 

For example, consider the source mass-multipole mo- 
ments written in the form 

h rn oc rjMr(t) l e-' lmv ^[l + 0(2)]. (5.10) 



r(l 



C. Memory in numerical simulations 

While post-Newtonian calculations can only reliably de- 
termine the memory during the inspiral phase, in a quasi- 



15 For example, by parameterizing the orbit equations in terms of the 
true anomaly, one can show that at leading-PN order an extreme 
hyperbolic orbit with eccentricity eo > 1 and impact parameter 
b has linear memory 



AZ 



(2) 



8?r r\M £ 



b 

(+) 



[5(e=F^ 



(-)n 



where the argument of periastron lies on the +x axis. Here the 
linear memory arises from different incoming and outgoing phase 



angles tp 
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± arccos(— 1/eo). See also Eqs. (15) of Ref. [32 



which are equivalent to this result. 
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circular binary most of the memory will accumulate dur- 
ing the merger and ringdown when GW emission is great- 
est [IB] (see also the discussion in Sec. MI CI above). Nu- 
merical relativity (NR) simulations would thus seem to 
be an obvious approach for calculating the full nonlinear 
memory from coalescing compact binaries. However the 
accurate computation of the GW memory from NR sim- 
ulations faces several challenges. 

Restricting our discussion to nonspinning, quasicircu- 
lar binaries orbiting in the x-y plane, numerical relativity 
simulations can most accurately calculate the dominant 
I = m = 2 mode of the waveform. However no memory 
effect is present in this mode. The nonlinear memory is 
only present in the m = modes which have not been ex- 
amined (to my knowledge) in any of the published papers 
on NR simulations of quasicircular, inspir ailing binaries. 
These m = modes tend to be much smaller in magnitude 
and depend more sensitively than other modes on the ini- 
tial conditions of the simulations. To see why this is so, we 
begin by considering the two ways in which gravitational 
radiation is typically calculated in NR simulations. The 
most widely used approach is to calculate the curvature 
scalar ^4 and decompose its value at large R into a sum 
over spin- weighted spherical- harmonic modes: 16 

00 / 

* 4 = h + -ih x =J2 J2 *Plm(t,R)-2Y lm {G,$). (5.12) 

1=2 m=-l 

Comparison with Eq. (|2.6j) shows that the modes of $4 
can be directly related to the modes of the metric pertur- 
bation: 

i>im = hi m . (5.13) 

Some NR simulations [H, [56|, [97], [IU also compute quan- 
tities which are more directly related to the h\ m metric 
perturbation modes [see Eq. (30) of Ref. [56] for the exact 
relation] . 

We can use PN waveforms to gain insight into the rel- 
ative magnitudes of the hi m and tpi m modes. The hi m 
modes have been calculated to 3PN order in Refs. [40. l84j ■ 
The m = pieces of hi m are computed to 3PN order in 
Sec. lIVI above. From these hi m modes it is straightforward 
to compute the ipi m modes. This is detailed in Appendix 
iBl Using these results, we plot in Fig. [2] the absolute val- 
ues of him and ipi m for some selected modes as a function 
of the PN parameter x. 17 Both plots show some of the 
largest oscillatory modes for rj = 1/4 — the (2,2), (4,4), 
(3, 2), and (4, 2) modes — as well as selected m = "mem- 
ory" modes: (2,0), (3,0), and (4,0). Notice that when 
the him modes are compared, the /120 mode has the sec- 
ond largest magnitude after the ft.22 mode. Note also the 



Note that different conventions exist in the literature for the re- 
lation between ^4, and h x - 

Similar plots from NR simulations that compare the evolution of 
different (I, m ^ 0) modes during the merger an d rin gdown can 
be found in Fig. 3 of Ref. jM, Fig. 11 of Ref. [Iq3, Fig. 2 of 
Ref. and Fig. 24 of Ref. |60|. 



PN scaling of these hi m modes: 

I&22I ~ \hm\ ~ M ~ 0(r]x), (5.14a) 

M ~ \hw\ ~ M ~ 0{nx 2 {l - 37 7 )], (5.14b) 

|/i 30 | ~ 0( V 2 x 7 / 2 ). (5.14c) 



Sill, even though the modes corresponding to the leading- 
order Christodoulou memory — /120 and /140 — enter at the 
same PN order as /122, they are significantly smaller in 
magnitude: | ^22 1 ~ 9.7|/i2o| ~ 850|/i4o| at x f» 1/5. This 
several orders-of-magnitudc difference between the dom- 
inant (2, 2) mode and the leading-order memory modes 
means that the memory effect will be difficult to resolve 
even for simulations that directly compute the modes of 
the metric perturbation. 

As the right plot of Fig. [5] shows, the situation 
is significantly worse for simulations that extract the 
gravitational-wave content from the 0; m modes. In this 
case the largest memory modes are many orders of mag- 
nitude smaller than the largest oscillatory modes. For 
comparison, the PN scaling of the ipi m modes is 



IV22I ~ 0(r)X 4 ), (5.15a) 

|-04 4 | ~ |032| ~ |-042| ~ 0[ V X 5 (1 - 377)], (5.15b) 

IV20I ~ hM ~ 0( V 3 x 9 ), (5.15c) 

|03o| - 0(rfx 23 / 2 ). (5.15d) 



In terms of the ipi m modes, the leading-order 
Christodoulou-memory modes are smaller than the 022 
mode by five PN orders. The relative mode magnitudes 
are |t0 2 2| « 69O|0 2O | « 4.5 x 1O 4 |0 4O | at x w 1/5. Signifi- 
cant improvements in the accuracy of NR simulations will 
be required to resolve these memory modes. Comparison 
with Fig. 2 of Ref. [6!| suggests that current simulations 
can resolve "0im modes that are about two orders of mag- 
nitude smaller than the dominant 1P22 mode. 

Numerical relativity simulations face further difficulties 
in calculating the memory. For those simulations that 
compute the tpi rn curvature perturbation modes, two time 
integrations are required to construct the hi m metric per- 
turbation modes. As discussed in detail in Ref. (69[ , choos- 
ing these integration constants to be zero causes an artifi- 
cial memory and a slope in the hi m modes (see their Figs. 4 
and 5). This artificial memory is completely nonphysical 
and appears for all (l,m) modes. It arises from the finite 
size of the initial separation and extraction radius. The 
initial burst of "junk" radiation could also contribute to 
this artificial memory. (These effects are likely respon- 
sible fo r the memory seen in the "Lazarus" simulations 
of Ref. [l0l| -) Reference [6!| suggests choosing one inte- 
gration constant so that the him modes have zero slope, 
and choosing the other constant so that the hi m modes 
have zero offset at late times. This procedure should pro- 
duce physically accurate metric waveforms for the m 7^ 
modes (which do not possess any physical memory). How- 
ever, this procedure will completely cancel any physical 
memory effect present in the m — modes. Since the ipio 
modes are so small anyway, this will not be an important 
issue until NR simulations become more accurate. When 
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FIG. 2: (color online). Comparison of memory and nonmemory modes. The left plot shows the absolute value of some of the 
largest hi m modes — as well as the largest m = modes — as a function of the post-Newtonian (PN) parameter x. The m 7^ 
modes are the solid (red) curves. From top to bottom they are /122, h.44, /132, and Kas,. These are computed from Eqs (9.3)-(9.4) 
of Ref. [40|. The remaining curves are the m = nonoscillatory modes. The long-dashed (blue) curves are the /120 (top) and 
hio modes [Eqs. (|4.3[l ] expanded to 3PN order. The short-dashed (navy) curve at the bottom is the /130 mode [the nonlinear, 
nonhereditary DC term given by Eq. (|5.8p ]. The right plot is similar except it shows the absolute value of the corresponding 
iplm modes (see equations in Appendix [B]). Note how the relative values of the oscillatory (m 7^ 0) and the memory (m = 0) 
modes change in the two plots. Although the memory is relatively large in the metric-perturbation modes, the memory modes 
are significantly suppressed relative to the other curvature-perturbation modes. Both plots are for equal-mass binaries (77 = 1 /4) 
and span the range x = 1/30 — 1/5. Recall that larger values of x correspond to smaller orbital separations and later times. Note 
that the modes with odd m vanish for equal-mass binaries. Several other numerically smaller modes would also appear on these 
plots, but are suppressed for clarity. 



that time arrives, one approach to computing those inte- 
gration constants is to simply match the hi m modes com- 
puted from the simulations with their corresponding PN 
expansions [Eqs. l|4.3jl ] at large orbital separations. 

Simulations that can compute the hi m modes directly 
do not have the problem of choosing two integration con- 
stants. They implicitly assume that the initial value of 
the asymptotic metric is zero at the start of the simu- 
lation. However, in addition to the problem of the /140 
and higher-order memory modes being small in magni- 
tude (and thus hard to resolve), these simulations (as 
well as those that compute ipi m ) must also deal with the 
memory's sensitivity to the past history of the source. 
To see how this sensitivity comes about, let us examine 
the leading- (Newtonian)-order piece of the largest mem- 
ory mode /i2o : 



"20 



(Tr) 



V2R 



Tr 



U. 



T 



(mcm)(l) 
20 



(t) dt, 



(5.16) 



where t/20 ' is given by the leading-order piece of 
Eq. (|3.13ap . and To is the starting time of the simulation. 
For a realistic, quasicircular binary that inspirals from a 
large initial separation, To — > —00. Using the Newtonian- 
order piece of Eq. (|3.19p , we can reexpress Eq. (|5.16[) as an 
integral from x (T ) « M/r (T ) to x(T R ) « M/r{T R ), 
yielding 
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Since this reduces to the memory from a realistic binary 
when ro — > 00 (or at least when r 3> M) , we can easily see 
that NR simulations of binaries at a finite initial separa- 
tion ro would underestimate the size of the memory. The 
resulting fractional error in the memory when the binary 
is at a smaller separation r is 
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For a simulation that starts with an initial separation of 
10M, the error in the memory at a harmonic coordinate 
separation of 5M is 50%. Thus, very large initial sep- 
arations (> 50M) are required to compute the memory 
accurately. (See also Sec. 4.2 of Ref. [35| for a related 
discussion.) As mentioned above, this issue could be cir- 
cumvented by using PN expressions for the memory as 
initial conditions at separations where PN theory is accu- 
rate. 

In addition to the problems mentioned above, other 
sources of error — such as gauge effects or a finite extrac- 
tion radius — can contaminate the numeric al w aveforms 
and possibly swamp a small memory signal [l02j |. 

While directly extracting the memory from the hi m or 
4>i m modes will be difficult, NR could be used to compute 
the memory indirectl y. T he effective-one-body (EOB) for- 
malism (see Damour 103] for references) provides a means 
to compute the h\ m modes semianalytically, calibrating 
their values to NR simulations. Reference 36] has used an 
EOB calculation of the h%% mode, combined with leading- 
order expressions for the energy flux, to compute the 
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memory for the merger and ringdown phases of binary 
black-hole coalescence. Another approach is to use a hy- 
brid PN/NR method to compute the memory. This would 
involve taking the numerical m ^ hi m modes (or the 
time integral of the ipi m modes) computed in an NR sim- 
ulation and substituting them directly into Eq. (|3.3[) . The 
resulting angular integrals could be evaluated analytically 
as was done in Sec. MI Al The final time integral must 
be performed numerically to compute the h^ e modes, 
with the PN expressions for the memory in Eqs. (4.3) used 
to compute the integration constant. This hybri d PN/NR 
procedure is currently being implemented [104] with the 
accurate merger waveforms described in Ref. [571 ]. 



D. Detecting memory with gravitational- wave 
interferometers 



While the gravitational- wave memory is an intrinsically 
interesting effect of general relativity, its scientific rele- 
vance relies on its potential for detection. The memory 
from a GW that passed through a region of space in the 
distant past (long before the start of the observation) is 
itself undetectable. In this case, the memory near a de- 
tector is simply a small, constant shift in the Minkowski 
metric, which is unobservable. For example, choosing a 
TT-coordinate system whose z axis is along the direction 
of the wave's propagation, the metric at late times is 



ds 



-dt A + {l + h + )dx z 



(1 



h + )dy 2 + 2h x dxdy + dz 2 , 
(5.19) 

where h+ and h x are constant in space and time near 
the detector. This metric has vanishing connection co- 
efficients (r^ 7 = 0) and can be transformed to a new 
coordinate system via 18 



T = t, 

Z = z, 

U = x- 



(5.20a) 
(5.20b) 



h x - h l l 2 J\-h\-h\ ) , (5.20c) 



V = xh 



1/2 y(h x h}_l 2 + ^i-hl-hl) 



I + Ha 



(5.20d) 



which has the flat metric 

ds 2 = -dT 2 + dU 2 + dV 2 + dZ 2 . 

What is actually observable is the buildup of the memory 
over some observation time as the GWs pass through the 
detector. 

The detectability of the memory was first studied by 
Thorne and Braginsky [!,[22| f° r unmodeled burst sources. 
Different types of detectors respond differently to the 
memory. For a memory with a characteristic rise time r, 



a resonant mass detector will have no significant response 
to the memory unless its resonant frequency is <C 1/r 
3]. A typical rise time for the Christodoulou me mory 
from a merging black- hole binary is r ~ (50— 100)M [36j |. 
Resonant mass detectors typically operate at frequencies 
> 700 Hz — far too high for them to be sensitive to the GW 
memory. An ideal laser interferometric detector whose 
test masses are completely unconstrained would experi- 
ence a permanent deformation after the GW has passed. 
In practice, a LIGO-type interferometer has servomech- 
anisms that keep the test mass positions fixed. So while 
LIGO is sensitive to the buildup of the memory, its design 
prevents it from physically "storing" the memory indef- 
initely. In contrast a detector like LISA is designed to 
be truly freely-falling and could (in principle) be perma- 
nently "deformed" by the passage of a GW with memory. 
However, we emphasize that this late-time deformation is 
not in itself directly observable except in comparison with 
some earlier nondeformed state. 

Properly computing the memory signal that a detector 
could measure requires a detailed knowledge of the mem- 
ory waveform throughout the inspiral and coalescence. 
This is particularly important for the memory (as op- 
posed to the oscillatory pieces of the GW) because most 
of the signal amplitude accumulates during the merger 
and ringdown phases. Kennefick [8(| has computed the 
signal-to-noise ratio (SNR) for the memory by truncating 
the Newtonian-order memory formula [Eq. (|4.6a[) ] at some 
critical radius r k /M > (QAnri) 2 ^ 5 w 4.8 for 77 = 1/4 [see 
his Eq. (21)] and assuming the memory remains constant 
at that value for later times. His approach has several im- 
portant limitations: (i) the Newtonian approximation is 
used in the strong field region where it is no longer valid; 
(ii) the final saturation value of the memory is set by a 
cutoff separation which is — to some extent — determined 
arbitrarily; and (iii) the rise and abrupt cutoff of the mem- 
ory signal does not necessarily mimic its true evolution. 
The details of how the memory rises and asymptotes to its 
final saturation value can significantly affect the Fourier 
transform of the memory signal and the resulting value 
for the SNR [cf. Fig. 2 of Ref. [H]. 

Rather than model the evolution of the memory sig- 
nal, the memory's SNR can be crudely estimated using 
the zero- frequency limit (ZFL) approximation fl4l . 1 105 " 



Il06l Il07|. This is essentially the approach followed by 
Thorne [221 ], who examined the Christodoulou memory 
from unmodeled GW bursts from merging BHs. 19 How- 
ever the ZFL is inadequate for computing the memory's 
SNR for some sources. As is also the case with Ken- 
nefick's analysis, the ZFL approach accurately computes 
the Fourier transform of the memory signal only in the 
regime where / < 1/r - 1/(70M) (36[. For two 10 6 M Q 
BHs, 1/r w 0.00145Hz; for two 10Af o BHs, 1/r « 145Hz. 
These frequencies lie near the peak sensitivity of space- 



If we require the coordinate transformation to be real- valued than 
we must impose the conditions /i+ > and + h 2 x < 1. 



Note that Thome's \'2"A normalization of the SNR should be cor- 
recte d by multiplying his Eq. (Al) by \/2 (see footnote 61 of 

Ref. fToai i. 
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based [8l| and terrestrial [4l| interferometric detectors, 
making the ZFL a poor tool for estimating the memory's 
SNR for these sources. Kennefick's [8(| estimate of the 
Fourier transform improves upon Thome's ZFL approach 
at low frequencies by incorporating a model for the mem- 
ory's rise during the early inspiral. But Kennefick's model 
also has large errors at frequencies ~ 1/t where much of 
the SNR is accumulated. Both approaches can potentially 
make order-of-magnitude errors in the SNR. An improved 
estimate of the memory's detectability needs to incorpo- 
rate a model for the memory's evolution during the merger 
and ringdown, as well as the inspiral. Such a model is de- 
veloped in Ref. [36J and compared with the approaches of 
Thorne [H and Kennefick ]m. 



VI. CONCLUSIONS 

The nonlinear (Christodoulou) memory is an interesting 
and potentially observable effect. Previous work has com- 
puted only the leading-order memory contributions to the 
gravitational- wave polarizations [H, [H, 53, • Here the 
post-Newtonian corrections to the memory for quasicircu- 
lar, inspiralling binaries have been computed. The main 
results are the equations listed in Sec. IIVI which gives the 
memory's contribution to the waveform's spin-weighted 
spherical-harmonic modes and the + polarization. These 
results are illustrated graphically in Fig. Q] and show that 
the post-Newtonian corrections to the memory are impor- 
tant in the late inspiral. The calculations presented here 
complete the memory to 3PN order. Other nonoscilla- 
tory effects on the waveform were also discussed. In par- 
ticular, a nonlinear, nonhereditary DC (nonoscillatory) 
term is present in the x polarization at 2.5PN order [3a |. 
while linear DC terms can enter the waveform at 5PN and 
higher orders. 

While most studies have focused on the memory in qua- 
sicircular orbits, the memory should also be calculated for 
other types of orbits. For hyperbolic orbits with small de- 
flection ang les, the nonlinear memory has been calculated 
in Ref. [32]. In future work the memory for inspiralling 
eccentric binaries [9l[ will also be computed. This is espe- 
cially important because circularized binaries have grow- 
ing eccentricity in the past, and the nonlinear memory is 
sensitive to the binary's past history. 

One of the main conclusions of this paper is that it will 
be difficult (but not impossible) to directly extract the 
memory from numerical relativity waveforms. Nonethe- 
less, numerical relativity will be critical in providing a 
full understanding of the memory. These simulations will 
eventually better determine the size of the memory in the 
late inspiral, and they will also be able to compute the 
evolution and saturation value of the memory during the 
merger and ringdown phases of coalescence [361 ] . The for- 
mulas presented here could be compared with the inspiral 
phase of those future simulations. They could also be 



used to set the memory's initial value at the start of a 
simulation. 

Detecting the memory in binary black-hole coalescences 
could serve as a test of general relativity. However, the 
nonoscillatory nature of the memory makes its detection 
difficult. Detections will only be likely for sources whose 
primary waves will have high signal-to- noise ratios, such 
as supermassive black-hole binaries in the LISA band. 
This paper argues that previous signal-to-noise ratio esti- 
mates of the memory may contain large errors. This issue 
is investigated further in Ref. (36|. 

Aside from the formal mathematical descriptions in 
Refs. [26l 1 109l | . the memory has only been explored in de- 
tail with post-Newtonian theory. It would be interesting 
to also explore the memory using second-order black-hole 
perturbation theory. Recent work along these lines has 
demonstrated a memory effect in the second-order quasi- 
normal modes of the Schwarzschild spacctime |110j - It 
would also be interesting to see how the nonlinear mem- 
ory manifests itself in alternative theories of gravity. 
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APPENDIX A: ANGULAR INTEGRAL OF THE 
TRIPLE PRODUCT OF SPIN- WEIGHTED 
SPHERICAL HARMONICS 

In this appendix we briefly derive a general formula for 
the angular integral of a product of three spin-weighted 
spherical-harmonic functions. Integrals of this form are 
needed when we evaluate the nonlinear memory contri- 
bution to the radiative mass-multipole moments Ui m [see 
Eq. (|3~3"]) ]. We begin by rewriting Eq. (|2~Ti]) as 



(©) = sA k i) sin - ] ( cos 



(Al) 



where 
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n(k)= (-IP [fe + mjWi rnWi + *j)Kh ~ *j)'] V2 r A2] 

9A 3> m, + - - Sj - k 3 )\{ Sj - mj + k 3 )\ ' 1 1 

Pj = 2kj + Sj — rrij, and k^/j) and fc/(j) are defined analogously to the limits in Eq. (|2.14[) . The index j = 1, 2, 3 serves 
only to label the three harmonics. The integrals we wish to evaluate can then be written as: 



~ (47r) 3 / 2 Jo J "' 



'» : *i " sin9 rf9 . 



(A3) 



The $ integral is simply 
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giCnu+ma+m,)* rf$ = 27^""^ . (A4) 



The 9 integral can be written as 

/ < Sl ^ 2S2 < S3 sin9de = 2 J2 9i(h)9 2 (k 2 )g 3 (k 3 ) sin- (cos- dQ , (A5) 



where a = 1 + (pi + p 2 + Ps)/^ and b = 1 + h + l 2 + I3 - 
(pi +P2 +p3)/2. The remaining O integral is expressible in 
terms of the Beta or Gamma functions and can be found 
in standard tables 



e 



sin ■ 



2a-l 



6 



cos ■ 



26-1 



r( a )r(6) 

This solution is easily extended to arbitrary products of 
the spin- weighted spherical harmonics. Integrals involv- 
ing the complex conjugates of one or more of the ~ s Y lm 
in the integrand can be computed by using s Y lm * = 
(— l) s+m ^ s Y l ~ m . While computer algebra programs such 
as Maple or Mathematica can symbolically evaluate the 
angular integrals in Eq. (|A3|) directly it is significantly 
faster for those programs to evaluate the result when ex- 
pressed as the above sum over Gamma functions. This is 
helpful when many such integrals need to be evaluated for 
large I values. 

We also note that the integral (|A3[) is proportional to 
the product of two 3-j symbols or two Clebsch-Gordon 
coefficients. This follows from the fact that the spin- 
weighted spherical harmonics are related to the Wigner 
matrices, which have a w ell-kn own triple-product inte- 
gral [cf. Eq.(4.6.2) of Ref. [H3, Eq. (4.62) of Ref. [IT3 |, 
Eq. (110.3) of Ref. or Eq. (1.42) of Ref. [lT|. How- 

ever, one must use extreme care when comparing nota- 
tion, definitions, and phase conventions among different 
authors. The above formulas provide an explicit result in 



terms of the Gamma function, which is uniformly defined. 
APPENDIX B: CALCULATION OF THE h lm AND 
i, lm MODES 

In this appendix we outline the computation of the ex- 
plicit PN expansions (for quasicircular orbits) of the ipi m 
modes defined in Eq. (|5.12p . These expressions are used 
to justify the results of Sec. IV CI and in generating the 
curves in Fig. [2j The %j)i m modes are directly related to 
the metric perturbation modes via ipi m = hi m . The /i; m 
modes are written as 



(Bl) 



where the Hi m are given up to 3PN order by Eqs. (9.4) 
of Ref. [!(| and Eqs. ()4.3j) above. The phase variable ^ is 
related to the orbital phase ip via [Eq. (8.8) of Ref. [401 ]] 



ip = <p- 





"1 ^ " 











(B2) 



where lnaio = i| — |te — fin 2 + | ln is related 

to the arbitrary constant tq appearing in the coordinate 
transformation in Eq. (|2.20p . The orbital phase can be 
expressed explicitly as a function of time, or more con- 
veniently as a function of x, to 3.5PN order [Eq. (5) of 
Ref. U]: 
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where we have substituted the values A = -1987/3080 and 9 = -11 831/9240 for the ambiguity parameters [9ll93l. l94j . 
For reference, we also list the full I = m = 2 mode [Eq. (9.4a) of Ref. [40j|]: 
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where a is +1 for the polarization conventions in Kidder [84J and —1 for the conventions in Blanchet et al. |40 |. 

To compute the time derivatives of the metric perturbation modes, we express hi m and ip entirely in terms of x using 
the above formulas. We then use Eq. (I3.19|) to express the time derivatives in terms of x derivatives: f = df/dt — if, 
where /' = df/dx for an arbitrary function f[x(t)]. A second-time derivative is then given by / = x(f"x + fx'), where 
x' = [right-hand side of Eq. (|3 . 19[) ] . Defining Hi m = xHi m , the first-time derivatives of the hi m modes are easily 
seen to be 
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These expressions are used in evaluating the energy flux [Eq. (|2.28p ] and the memory contribution to the radiative-mass 
multipoles [Eq. (|3.3[) ]. For reference we list only the value for the I = m = 2 mode: 
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The values of the hio modes can be inferred from Eqs. (|2.7|) and (|3.13p above. 
Taking another time derivative lets us write the ipi m modes as 
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i>l m = k m = 8J^ 1 ^-xe- i "'' 



(-imii')H lm + H' lm +x Hl' m -2imip'H' lm -m(iip" + mi\)' 2 )U lm }, 



V 5 MR Vim 



(B7) 



We list here the I = m = 2 and the m — modes 



2(1 



Note that the expression for ip22 agrees with the 2.5PN order 
result in Eqs. (3.7a) and (Ala) of Ref. fsotl when their parameter 

■uj = —24. Note also that the 3PN order expressions for -H22, ^22, 
and i/>22 are almost identical except for the imaginary part of the 
2.5PN order term. 
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The remaining tpi m modes for m ^ can be easily computed by applying Eq. (|B7[) to Eqs. (9.4) of Ref. [4(| (see also 
Appendix A of Ref. (69] ) ■ These modes (including some not listed here) are used in the construction of Fig. [2] 
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